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The negative sign problem in quantum Monte Carlo (QMC) simulations of cluster impurity prob¬ 
lems is the major bottleneck in cluster dynamical mean field calculations. In this paper we sys¬ 
tematically investigate the dependence of the sign problem on the single-particle basis. We explore 
both the hybridization-expansion and the interaction-expansion variants of continuous-time QMC 
for three-site and four-site impurity models with baths that are diagonal in the orbital degrees of 
freedom. We find that the sign problem in these models can be substantially reduced by using 
a non-trivial single-particle basis. Such bases can be generated by diagonalizing a subset of the 
intracluster hoppings. 

PACS numbers: 71.10.Fd 


I. INTRODUCTION 


Quantum Monte Carlo (QMC) methods are a powerful 
tool for studying the properties of quantum many-body 
systems. They are based on a mapping of the quan¬ 
tum system to an auxiliary classical system, which is 
then sampled stochastically. A fundamental limitation of 
fermionic QMC algorithms is the so-called negative sign 
problem,® which appears when configurations with neg¬ 
ative weights appear due to fermionic statistics. Apart 
from models with special symmetries,®® the sign prob¬ 
lem prevents access to the low-temperature bulk proper¬ 
ties of many-body fermionic systems since the efficiency 
of the MC sampling decreases exponentially with system 
size and inverse temperature as detailed below. 

In a QMC algorithm, we express the partition function 
of the quantum system as the sum of contributions of 
classical configurations {c}, which define a configuration 
space Q, 

Z = Tr[e~^ H ] = Y, w - (!) 

cefi 


The expectation value of an observable O is then given 
by 


<o) = 


Tr [Oe-P H ] 
Tr[e~P H ] 


= 1 X °c w c 

cen 


( 2 ) 


where O c is the value of the observable associated with 
c. When w c > 0 (Vc), ( O ) can be efficiently evaluated 
by a summation over N c configurations sampled from f 1 
according to the weight {w c }, 

(O) ~ (0) MC = . (3) 

The standard deviation 5 ( O ) of this MC estimate scales 
as 1 /y/N~ c when N c is larger than the autocorrelation 
time of the MC dynamics. 


The MC sampling suffers from a negative sign problem 
when w c < 0 for some configurations. In this case, one 
cannot interpret w c /Z as the probability for the configu¬ 
ration c. Using |w c | for the weight instead, one can still 
perform importance sampling as follows: 

(Ay. = E c sig n (c) O(c) (sign Q) mc 
E c sign(c) (sign) MC 

However, in simulations with a sign problem, both 
(sign 0) MC and (sign) MC decay exponentially with in¬ 
creasing system size and inverse temperature. This 
makes the MC sampling exponentially inefficient. To 
keep the error on (O) constant, we need to increase the 
number of Monte Carlo steps (at least) as 1/ (sign)J^ c . 

In this paper we investigate the fermionic sign prob¬ 
lem for so-called quantum impurity models. 10 These de¬ 
scribe a small number of interacting orbitals hybridized 
with a noninteracting bath. The simplest of these mod¬ 
els is the single-orbital Anderson impurity model which 
was originally propose d to describe a magnetic impurity 
embedded in a metal. 11 12 Subsequently, it was found 
that in infinite dimensions, the Hubbard model can be 
exactly mapped onto a single-site impurity model with 
a bath which is self-consistently determined ! 13 ! 14 ! In fi¬ 
nite dimensions, the analogous procedure leads to the 
dynamical mean-field theory (DMFT) approximation for 
correlated lattice models.® 15 DMFT can be extended to 
incorporate short-range correlation^® by mapping onto 
cluster-type imp urity systems, as done in schemes such 
as cluster DMFTp®33and the dynamical cluster approx¬ 
imation.®^ The success of DMFT created a demand for 
powerful and versatile impurity solvers and led to the 
development of continuous-time Monte Carlo impurity 
solvers. 

There are two complementary types of continuous-time 
quantum Monte Carlo algorithms for solving an impu¬ 
rity problem, which are both based on a stochastic sam¬ 
pling of a perturbation expansion: the interaction ex- 
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pansion method (CT-INT)J^ the auxiliary-field Monte 
Carlo meth od and the hybridization expansion method 
(CT-HYB) ! 22 1 2 3 1 These methods do not have a sign prob¬ 
lem at an y filling in the case of the single-orbital impu¬ 
rity model. 22 24 25 -l However, they generally suffer from a 
negative sign problem when applied to a multi-orbital 
or cluster impurity model with an internal structure in 
which electrons can be exchanged and the sign problem 
becomes worse as the system becomes larger. 

Since the sign problem is not gauge invariant and rep¬ 
resentation dependent, the severity of the sign problem 
depends on the single-particle basis used to represent the 
impurity. For CT-HYB, it has been empirically known 
that the sign problem is sometimes improved by using 
the single-particle basis that diagonalizes the intraclus¬ 
ter single-particle Hamiltonian of the impurity ! 26 ^ 27 How¬ 
ever, to the best of our knowledge, no systematic effort 
has yet been made to clarify how the sign problem de¬ 
pends on the single-particle basis for CT-HYB. On the 
other hand, for CT-INT, the site basis, with local interac¬ 
tions identical to those of the lattice model, is commonly 
used in solving Hubbard-like cluster impurity models. 
There, one faces a severe sign problem particularly in 
systems away from half filling and in geometrically frus¬ 
trated systems. However, there has been no previous 
effort to test alternative single-particle bases to reduce 
the sign problem. 

In this work, we systematically investigate the single- 
particle-basis dependence of the sign problem for typical 
cluster impurity models within CT-HYB and CT-INT. 
We focus on a trimer model and a tetramer impurity 
model with diagonal baths. The main conclusion of this 
study is that the sign problem can be dramatically im¬ 
proved by choosing a single particle basis which diagonal¬ 
izes some part of the intracluster single-particle Hamil¬ 
tonian. 

The rest of this paper is organized as follows. Sec¬ 
tion II is devoted to the results obtained for CT-HYB. 
Section II A describes the formalism and the implemen¬ 
tation of the basis rotation. We present the results for 
the trimer and the tetramer models in Secs. II B and II C, 
respectively. Section III discussed the basis dependence 
for CT-INT. The implementation is described in Sec. Ill 
A. Then, we show the results of the trimer and tetramer 
models in Secs. Ill B and III C, respectively. A summary 
and discussion are provided in Sec. IV. 


II. HYBRIDIZATION-EXPANSION 
A. Implementation 


by 


^_, rP 

S = S loc + ^2 dTdT'A afe (T' - t)c4(t')q,(t), (5) 

ab J ° 

where A is the hybridization function [A a f ,(t) = AjJ a (r)] 
and a, b are combined spin and site indices. We denote 
the inverse temperature by /3. If S\ oc contains only in¬ 
stantaneous terms up to two-body interactions, the cor¬ 
responding Hamiltonian reads 

%loc — ^ ' tabCa^b Y ^ ) ^abcdd^C^CcCd- (6) 

ab abed 

We may define a transformed single-particle basis 

d a =J2 U baCb, (7) 

b 

4 = ( 8 ) 

b 

with U a b being a unitary matrix. In this basis the second 
term in Eq. 0 reads 

r P 

>Shyb = E/ drdr,A o^-( r, “ T ) d l( T ')db(T), (9) 

ab 

with 

A ab (r) = E(C7 t )acA cd (r)C/ d6 . (10) 

cd 


Now, we expand the partition function Z as 


Z - Z hath E n|2 E E 

n=0 ’ 

rP rP 

/ dridr[ • • • / dry dr' 

Jo Jo 

Tr loc e~ / 3 Wloo Td an (T n )dE T n) ’ ’ ’ d ai (n)d^ (r[) 

( 11 ) 


xdet M 1 


where Zbath is the partition function of the bath. The 
matrix elements of (M )j iJ are given by the hybridiza¬ 
tion function A Q ,' iQ .(r J ' — tj). 

The partition function in Eq. (11) can be evaluated 
by importance sampling of configurations of annihilation 
and creation operators on the imaginary-time axis. The 
weights are given by 


w(d ai (n),--- , d arl (T n ); da^r'O, ■ ■ ■ ,d a '(' r' n )) 
e~ PHloc Td an (T n )dL «)•• 


dr 


-Tr 


d ai (Ti)d\' (r[) det M \ 


( 12 ) 


In this section,we review the hybridization expansion 
continuous-time quantum Monte Carl o algo rithm (CT- 
HYB) for quantum impurity problems ! 22 * 2 ^ Tracing out 
the bath degrees of freedom, the effective action is given 


The local trace is evaluated using the matr ix f ormalism 
where the operators e~ rWloc , d a , in Eq. (12) are rep¬ 
resented in the eigenbasis of 7di oc - Note that we do not 
have to transform Hioc to the new single-particle basis. 
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B. Trimer model 


1. Set-up 


In this section, we study a trimer impurity model as a 
minimal system which exhibits a negative sign problem. 
As shown in Fig. [Tj each site is connected to a bath with 
a semicircular density of states of width 4. The Hamil¬ 
tonian is given by 

3 3 3 

H = - E ^ E ~ /^E^t 

*=i *=i 

+A EE ( c lcr 0 afco- + a akcr C i•&) 

i—ot kcr 

+EE e k a afccr a a/c<T j (13) 

a: kcr 


with pL = U /2 and A = 1. Throughout this section, we 
measure energy in units of A. We set the onsite Coulomb 
repulsion to U = 5 throughout this section. The intr¬ 
acluster hopping matrix elements are given by f 12 = if , 
*13 = *, *23 = t' (kj = tji ) [see Fig. JTfa)]. a is the in¬ 
dex of the bath and the different bathlevels are labeled 
by k. The distribution of their energy levels ej~ is given 
by Ek S ( e ~ e fc) = f+V^l — e 2 , that is, the hybridization 
function A is given by 


AMr) = f E 1 

TI— — 00 


- 


1-2 


X 00 

= -J- E^ w ™ - ^+4) Sin(w„r) (0 < r < /3), 

P n—0 

(14) 

where co n = (2n + l)7r//? is the Matsubara frequency.^ 
Note that the hybridization function is invariant with 
respect to single-particle basis rotations. 

Since the hybridization function is diagonal, the nega¬ 
tive sign problem arises from the competition between t 
and U. The system exhibits no sign problem in the two 
limiting cases t = 0 and U = 0. In these two cases, the 
impurity problem is diagonalized by choosing U to be 
the identity matrix or the matrix made of eigenvectors 
of 'Himp, respectively. We call the former basis the site 
basis and the latter basis the diagonal basis. 

To find the optimal basis in terms of the average sign, 
we perform a brute-force search in the parameter space of 
U, restricting ourselves to (real) orthogonal matrices £ 
SO{ 3) which can be parameterized by three Euler angles. 
The average MC sign is computed on a uniform grid of 
20 x 20 x 20 for the Euler angles. Throughout this section, 
we set U = 5 and vary t and t'. 


2. Results for t = t' 

First, we discuss the results for t = if . Figure [2] shows 
the average sign computed for the site basis and the diag¬ 


onal basis as a function of t for U = 5 (/? = 15, 25, 50). 
We also show the highest average sign found in the brute- 
force search for several parameters. At j3 = 15, the aver¬ 
age signs for the site and diagonal bases cross at t ~ 0.7. 
Although the sign problem is severe for the site and di¬ 
agonal bases at * ~ 0.7, the optimal basis found by the 
brute-force search has a considerably higher average sign 
with weaker f3 dependence. We found that the optimal 
basis corresponds to the transformation matrix 


U = 


i 



v*2 

0 



(15) 


and therefore consists of bonding and anti-bonding or¬ 
bitals on one of the three edges (1 and 2) and a localized 
orbital on the remaining site (3). We call this basis the 
dimer+monoiner basis. (There are three equivalent op¬ 
timal bases for t = t'.) 

Figure [2] also shows the average sign for the 
dimer+monomer basis computed for 0 < t < 1. This ba¬ 
sis has a higher average sign than the site and diagonal 
bases in this entire parameter region. Another interesting 
observation is that the dimer+monomer basis is negative- 
sign-free in the two limiting cases: t U and t = 0. 
This may be because the system has no fermionic loop in 
the dimer+monomer basis for these two limiting cases. 
As illustrated in Fig. 0M. there is no hopping between 
the bonding and anti-bonding orbitals by construction. 
Instead, the basis rotation generates off-site Coulomb in¬ 
teractions between them (including non-density-density 
terms). The explicit form is given in Appendix [D] When 
t = 0 and/or U = 0, the system contains no fermionic 
loop as is clearly seen in Fig. 0b). This gives an impor¬ 
tant insight, namely that the optimal basis diagonalizes 
some subset of the intracluster hoppings. This allows to 
interpolate between the site and diagonal bases which are 
optimal in the limits * = 0 and U = 0. 

Figure[2](b) shows the /3 dependence of the average sign 
computed for the site basis, the diagonal basis, and the 
dimer+monomer basis at * = 0.6. These data fit the 
exponential form 

(sign) oc exp (■--E) , (16) 

V Psign J 

with /3 s j gn = 28 + 1, 40 ± 2, 192 ± 13, respectively. For 
the dimer+monomer basis, the average sign (sign) decays 
about five times more slowly with respect to /3 than for 
the diagonal basis. 

In Fig.|3j we compare the imaginary-time Green’s func¬ 
tion G(t ) computed for t = 0.5 and (3 = 50 using the site 
basis and the dimer+monomer basis and the same num¬ 
ber of MC steps. We plot the on-site (diagonal) element 
of G(t) in the site basis. Both results agree within sta¬ 
tistical errors (indicated by the noise), but the errors of 
the dimer+monomer basis simulation are much smaller 
due to a larger average sign. 

Figure [4] shows the expansion-order-resolved average 
sign (sign) ra computed at t = t' = 0.6 and /3 = 
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(a) Bath 




(c) N l =3, No. 1 


50 using the site basis, the diagonal basis, and the 
dimer+monomer basis. The distribution function of the 
expansion order has a single peak around n = 85 indepen¬ 
dently of the choice of the single-particle basis. On the 
other hand, the average sign (sign) n decreases monoton- 
ically for n < 85 for these three bases, becoming almost 
independent of n around the peak of the distribution 
function. The superiority of the dimer+monomer basis is 
already discernible at low expansion orders n — 20. Note 
that the distribution function P(n) generally depends on 
the single-particle basis used to represent the impurity. 
In Fig. [4j however, P{n) seems to be almost identical 
for these three bases. As detailed in Appendix [C| P(n) 
is representation independent when (sign}„ does not de¬ 
pend on n even if (sign)„ itself depends on the single¬ 
particle basis. Our observation for P(n ) is consistent 
with the small n dependence of (sign) n around n ~ 85. 



N l =0, No. 1 


• • 

FIG. 1. (Color online) (a) Trimer impurity model with 
orbital-diagonal bath. Each site is connected to a bath with 
a semicircular density of states of width 4i. (b) Schematic 
representation of the hoppings and the Coulomb interactions 
in the dimer+monomer basis (t = t'). (c) Inequivalent real¬ 
izations of links for the trimer impurity problem (t ^ t'). Nl 
is the number of links in a graph. 


3. Results for t 


We next discuss the results for t ^ t '. The result for 
t = t' suggested a guiding principle for generating good 
single-particle bases: diagonalizing a subset of the local 
hoppings. We thus introduce the concept of partial diag- 
onalization , i.e. we define a new basis as the eigenbasis of 
some selected matrix elements of the intracluster single¬ 
particle Hamiltonian. In practice, the trimer impurity 
model has three edges (links) with local hoppings. We 
deactivate some of them for the purpose of partially di¬ 
agonalizing the intracluster single-particle Hamiltonian. 
Figure nTc) shows all symmetrically inequivalent graphs 
for t 7 ^ v . Nl is the number of links remaining in a graph. 
The graphs for Nl = 0 and 3 give the site basis and the 
diagonal basis, respectively. Two dimer+monomer-type 
bases are generated from Nl, = 1 graphs. The eigenbases 
are generated by diagonalizing the intracluster single¬ 
particle Hamiltonians of the graphs numerically. Since 
some of the graphs have degenerate eigenvalues, we spec¬ 
ify the unitary matrices used for the present study in 
Sec. |A] 

We computed the average sign for these bases at t'/t = 
0.5 and 2.0 (/3 = 50). Figures [5ja) and J5|b) show the 
highest average sign among the graphs with the same 
Nl for t'/t = 0.5 and 2.0. For the site and diagonal 
bases, the sign problem is severe around t = 0.4 and 1.0, 
respectively. In these regions, the highest average sign is 
produced by Nl = 1 (No. 2), where a dimer is placed on 
the bond connecting sites 1 and 3. This basis is always 
superior to Nl = 1 (No. 1). 
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(a) 







FIG. 4. (Color online) Distribution function of the expan¬ 
sion order (n) and expansion-order-resolved average sign com¬ 
puted at t = t' = 0.6 and fi = 50. 


C. Tetramer impurity model 

1. Set-up 



The tetramer impurity model is illustrated in Fig. [6] 
Its Hamiltonian reads 

4 4 4 


^ = ~ E tijcl^Cja + UYl n it n 4 - /i E hi 

i^j i—1 i =1 

+ ^ E/ ^^( c ia a aka + a akcr C i<r) 

i—cx ka 

+EE e k a ]xkcr ao ‘kcn (17) 

cl ka 


FIG. 2. (Color online) Average sign computed for the trimer 
model with t'/t = 1 and U = 5 (ft = 15, 25, 50). (a) t depen¬ 
dence of the average sign. We show data for the site basis, 
the diagonal basis, the optimal basis found in the brute-force 
search, and the dimer+monomer basis, (b) /3 dependence of 
the average sign computed for the site basis, the diagonal ba¬ 
sis, and the dimer+monomer basis at t = 0.6. The dotted 
lines denote a fit by the exponential form (161. 



FIG. 3. (Color online) On-site imaginary-time Green’s func¬ 
tion G(r) computed for the trimer model with t = t' = 0.5 
and U = 5 (P = 50). We symmetrize the data using the 
three-fold rotational symmetry. 


where /r = U/2 and A = 1. Throughout this section, we 
measure energy in units of A. As in the case of the trimer 
impurity model, each site is connected to a bath with a 
semicircular density of states of width 4. That is, the 
distribution of their energy levels e*, is given by <5(e — 
£fc) = 2 + \/4 — e 2 - The definition of the matrix elements 
{tij} is shown in Fig. When t = t', this model is 

equivalent to a tetrahedron impurity model with cubic 
symmetry. The ratio t'/t (< 1) controls the strength of 
the geometrical frustration. Throughout this section, we 
set U = 5 and vary t and t'. 


2. Results for t = t 1 

Let us start by considering the results for t = t' . Fig¬ 
ure pla) shows all symmetrically inequivalent graphs for 
t = i. We obtained unitary matrices by numerically di¬ 
agonalizing the corresponding intracluster single-particle 
Hamiltonians.^ Some graphs have degenerate eigenval- 
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FIG. 5. (Color online) Average sign computed for the trimer model: (a) t'/t = 2 and (b) 0.5 at /? = 50. We show the maximum 
value of the average sign among single-particle bases generated by graphs with the same Nl■ We also plot the highest value of 
the average sign found in a brute force search. We show the average sign computed for each basis in Fig. [Tl]of the Appendix. 


ues and hence the unitary matrix is not uniquely deter¬ 
mined. For example, we obtained 


U = 


/ -1/2 
-1/2 
- 1/2 
V-V2 


3 

2 V3 



2\/3 


0 


1 



0 \ 

1 



1 4- 1 

2 ~ 2\/3 / 


(18) 


with the eigenvalues (—3,1,1,1) for Nl = 6 (No. 1). 
We provide all these unitary matrices in Appendix [A] In 
addition, we test the A\ (T 2 ) representations of the Tj 
point group 


U = 


1 

2 


11 1 1 \ 
1-1-1 1 I 

1 1 —1 —1 I ’ 
1-1 1 -1 / 


(19) 


which diagonalizes the intracluster single-particle Hamil¬ 
tonian of the tetramer as Nl = 6 (No. 1) does. We name 
this basis A\ (T 2 ). Nl = 0 corresponds to the site basis 
with U being the identity matrix. 

Figure [7)(a) presents the average sign computed at 
fj = 50. We show only the maximum values of the aver¬ 
age sign among single-particle bases generated by graphs 
with the same Nl- Refer to Fig. [12] in the Appendix for 
the data for each graph. The system exhibits a severe 
sign problem for the site basis [Nl = 0 (No. 1)], Nl = 6 
(No. 1), and Ai (T 2 ) at f ~ 0.4. A notable point is 
that Ai (T 2 ) is superior to Nl = 6 (No. 1) in the entire 
range of t. The only difference between these two bases 
is in the structure of the Coulomb matrix. At i ~ 0.4, 
the highest average sign is given by one of our nontrivial 
bases, Nl = 2 (No. 2) which consists of two independent 


dimers. Its unitary matrix reads 


1 V2 

0 

0 

n/2 \ 

0 

i 

V2 

1 

G2 

0 

1 

V2 

0 

0 

1 

\/2 

V 0 

1 

v/2 

1 

\/2 

0 / 


( 20 ) 


and the eigenvalues of the intracluster single-particle 
Hamiltonian of the graph are (—1, —1,1,1), respectively. 

On the other hand, the Nl = 4 (No. 1) basis has a 
rather high average sign comparable to A\ (T 2 ) overall. 
While this basis is numerically constructed to diagonalize 
the intracluster single-particle Hamiltonian correspond¬ 
ing to Nl = 4 (No. 1), it is found to also diagonalize that 
of Nl = 6, consisting of a different linear combination of 
the three-fold eigenvectors from Nl = 6 (No. 1) and Ai 
(T 2 ). Although one may be able to perform a brute-force 
search for the best linear combination, such a study is 
left for the future because of its high computational cost. 

Figure[7j'b) shows the f3 dependence of the average sign 
computed for the site basis, A\ (T 2 ), Nl = 2 (No. 2) at 
t = 0.4. The data fit the exponential form (16) with 
Aiign = 18.8 ±1, 48 ± 4, 67 ± 6, respectively. The average 
sign decays most slowly for Nl = 2 (No. 2) among these 
three. 


3. Results for t < t' 


For t < t ', there are more symmetrically inequivalent 
graphs, as shown in Fig. [6))b) , due to the absence of the 
cubic symmetry. The Nl = 6 (No. 1) basis diagonal¬ 
izes the local hoppings of the N L = 6 graph with the 
eigenvalues of (—2.5,0.5, 0.5,1.5). Some bases obtained 
for other graphs consist of different linear combinations 
of these doubly degenerate eigenvectors. For example, 
for the loop-type graph N L = 4 (No. 2), we obtained the 
same transformation matrix as that of A\ (T 2 ) [Eq. (19)]. 












7 
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FIG. 6. (Color online) Tetramer impurity model (inset) and symmetrically inequivalent graphs for t = t' [(a)] and t [(b)]. 


This transformation matrix also diagonalizes the intra¬ 
cluster single-particle Hamiltonian of N L = 6. 


We illustrate the t dependence of the average sign com¬ 
puted for if /t = 0.5 in Fig. [7](c). We show only the max¬ 
imum values among graphs with the same Nl- Refer to 
Fig. [T2] in the Appendix for the data for each graph. The 
trend is similar to t'/t = 1. For larger t, the highest av¬ 
erage sign is given by the loop-type graph Nl = 4 (No. 
2 ) ESI while the two-dimer-type graph Nl = 2 (No. 1) is 
optimal for smaller t. We plot the f3 dependence of the 
average sign at t = 0.4 in Fig. m- The data fit the 
exponential form (16) with /3 s ; gn = 21 ± 1, 44 ± 2, 78 ± 3 
for the site basis, Nl = 4 (No. 2), and N L = 2 (No. 1), 
respectively. 


III. INTERACTION-EXPANSION 

A. Implementation 


Motivated by the results for CT-HYB, we study the 
single-particle basis dependence of the sign problem for 
interaction-expansion continuous-time Monte Carlo (CT- 
INT) where we expand the partition function in powers 
of the interaction terms. We refer to Ref. [5T1 for a review 
of CT-INT. Even if the original model contains only on¬ 
site repulsion, a single-particle basis rotation may gener¬ 
ate non-density-density off-site interactions. The model 
presented in the new basis can be regarded as a multi- 




































(c) Tetramer model i!/t = 0.5 





FIG. 7. (Color online) Average sign computed for the tetramer model: (a)/(b) t'/t = 1 and (c)/(d) t'/t = 0.5. In (a) and (c), 
we show the maximum value of the average sign among single-particle bases generated by graphs with the same Nl computed 
at (5 = 50. In (b) and (d), we show the /3 dependence of the average sign computed at t = 0.4. The dotted lines are fits by the 
exponential (161. 


orbital impurity problem. However, for CT-INT, we have 
to expand in these interaction explicitly in contrast to 
CT-HYB. Therefore, in this paper, we restrict ourselves 
to single-particle bases for which non-local interactions 
are of “Slater-Kanamori” type, i.e., inter-orbital repul¬ 
sion, intra-orbital repulsion, exchange coupling, Hund 
coupling, and pair hopping. This can be done by re¬ 
stricting ourselves to graphs consisting of independent 
dimer (s) such as the dimer-|-monomer basis for the trimer 
and Nl = 2 (No. 2) for the tetramer shown in Fig. |6j(a) . 
For these bases, the Slater-Kanamori interaction with 
an unusual parameterization acts between each pair of 
a bonding orbital and an anti-bonding orbital as 

2 } ' ^ ' Ug^a'p , d\ ia c( i Q a ,dpi ( j'd a i a , ( 21 ) 

afta.'P' gg' 

with the intra-orbital repulsion U aoiaa = U / 2, the inter¬ 
orbital repulsion U a p a p = U' = U/2, the exchange in¬ 
teraction U a pp a = Jh = U/2, and the pair hopping 
Uaa /30 = = U/2. a and /3 are the bonding/anti¬ 

bonding orbitals (a ^ f3). cr and a' are spin indices. The 
derivation is given in Appendix [D] We treat these mod¬ 
els by an e fficien t algorithm recently proposed by some of 
the authoriplM) with the further optimization described 


in Appendix [E| While this transformation increases the 
number of interaction types, it does not increase the av¬ 
erage perturbation order since the strength of each in¬ 
teraction decreases. The Green function as a function 
of imaginary time is represented using 100 Legendre or¬ 
thogonal polynomials.^ 


B. Trimer 


We first study the trimer impurity model (131 for 
t = t' = 1. In Fig. §a), we show the U dependence 
of the average sign computed using the site basis and 
the dimer-l-monomer basis at (3 = 50. For the site basis, 
we face a severe sign problem for U > 7. On the other 
hand, the sign problem is substantially reduced for the 
dimer+monomer basis. The average sign is about 0.95 
even at U = 10 for this basis. 

Figure |8](b) compares the /3 dependence of the aver¬ 
age sign for these two bases at U = 8. The average 
sign vanishes exponentially with increasing /3. Fitting 
the data by Eq. (16), we obtained /3 s i gn = 59.5 ± 0.3 and 
1634 ± 17 for the site basis and the dimer+monomer ba¬ 
sis, respectively. This indicates that the dimer+monomer 
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FIG. 8. (Color online) Comparison of the average sign 
computed for the trimer model using the site basis and the 
dimer+monomer basis by CT-INT. (a) U dependence of the 
average sign at /? = 50. (b) /3 dependence of the average sign 
at U = 8. 




FIG. 9. (Color online) Comparison of the average sign com¬ 
puted for the tetramer using the site basis and the Nl = 2 
(No. 2) basis by CT-INT. (a) U dependence of the average 
sign at /3 = 50. (b) /3 dependence of the average sign at U = 8. 


basis allows to reach about 28 times lower temperatures 
compared to the site basis. 


C. Tetramer 


For the tetramer (t = t' = 1), we test Nl = 2 (No. 2) 
given in Eq. (201, which consists of two independent pairs 
of bonding and anti-bonding orbitals. In Fig. |9j)a) , we 
compare the U dependence of the average sign computed 
for the site basis and Nj j = 2 (No. 2) at /? = 50. For 
the site basis, one can see a severe sign problem for U > 
6. Note that the band width of the localized tetramer, 
which is not connected to the bath, is 6. On the other 
hand, for Nl = 2 (No. 2), the average sign decreases 
substantially more slowly. We plot the P dependence of 
the average sign at U = 8 in Fig. |9])b). The average 
sign vanishes following the exponential form (16) with 
p .sign = 20.58 ±0.07 and 154.2 ±0.9 for the site basis and 
Nl = 2 (No. 2), respectively. 

To demonstrate the advantage of Nl = 2 (No. 2), we 
compare the self-energy E(iw rl ) computed by using the 
site basis and the N L = 2 (No. 2) basis at U = 8 and 
P = 50. We used the same MC time for the site basis 
and Nl = 2 (No. 2) to make a fair comparison. Refer to 
Appendix |E| for technical details. The data obtained by 
the Nl = 2 (No. 2) are substantially less noisy for both 
of the site-diagonal and site-off-diagonal elements. 


IV. SUMMARY AND DISCUSSION 


In this paper, we showed that one can substantially 
reduce the negative sign problem in QMC simulations 
of quantum impurity models by using single-particle 
bases obtained by diagonalizing a part of the intracluster 
single-particle Hamiltonian. First, we investigated the 
trimer and the tetramer using CT-HYB. We found that 
optimal bases can be generated by diagonalizing subsets 
of the intracluster single-particle Hamiltonian in the im¬ 
purity. In particular, we revealed that single-particle 
bases consisting of bonding and anti-bonding orbitals 
have a high average sign when the kinetic energy and the 
onsite repulsion compete. Furthermore, we tested these 
bases in the framework of CT-INT for the trimer and 
the tetramer models. We showed that the sign problem 
is substantially suppressed with these bases, especially in 
the strongly correlated regime. 

In the present study, we restricted our consideration to 
orbital-diagonal baths to focus on the negative sign prob¬ 
lem arising from the local part of the impurity. Our new 
bases may be useful when one applies cluster extensions 
of DMFT to Hubbard models on frustrated lattices such 
as a Kagome lattic^SHdoJ an d a pyrochlore lattice.^ On 
the other hand, for CT-HYB, off-diagonal elements of the 
hybridization function also give rise to a sign problem. It 
remains to be clarified how to choose the optimal single¬ 
particle bases for impurity models with large off-diagonal 
elements. This will be practically important for investi- 
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lUJn 


ismP Although we employed the matrix formalism in 
the present study, the single-particle basis transforma¬ 
tion applies to the Krylov formalism as well. The Krylov 
algorithm may be more efficient than the matrix formal¬ 
ism when treating more than five orbitals. 

For CT-INT, we restricted our consideration to 
single-particle bases consisting of independent pairs of 
bonding/anti-bonding orbitals. The sign problem may 
be further reduced for more general single-particle bases 
which mix more than two sites. Such basis transforma¬ 
tions may generate general two-body interaction terms 
Vijkid\ a d' ja ,d ka 'd,i' T with i ^ j ± k ^ l (i, j, k, and l are 
orbitals in the rotated basis) such as correlated hoppings. 
One can still treat these general two-body interactions 
within CT-INT P 

Finding the best single-particle basis is essentially a 
classical optimization problem. Unitary matrices U can 
be parameterized as U = e lff , where H is a Hermitian 
matrixP We may be able to solve this optimization prob¬ 
lem numerically, e.g., by using simulated annealing ! 45 ! 46 ! 
Alternatively, constructions inspired from genetic algo¬ 
rithms might be useful: in such a “genetic” scheme, basis 
sets for smaller clusters might be used as “genes” from 
which trial basis sets for larger clusters are constructed. 
Mutations would then correspond to adding or removing 
links in the diagonalization process, and selection would 
be based on the average sign as the fitness function. In 
this way, one might hope to extend the present insights 
to clusters of much larger size than in the present study. 
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Appendix A: Single-particle rotation matrix 

We provide the single-particle rotation matrices used 
in the present study as text files. 


1. Trimer model 

The matrices for t'/t = 0.5 and t'/t = 2 are 
stored in “rotmat-trimer-tratioO.5.txt” and “rotmat- 
trimer-tratio2.0.txt”, respectively. The first, second, 
third columns are the indices i, j, and the entry £7^, 
respectively. The graphs are numbered from Nl = 3 in 
decreasing order and from No. 1 in increasing order as 
“Graph 1, 2, •••”. 


2. Tetramer model 

The matrices for t'/t = 1.0 and t'/t = 0.5 are 
stored in “rotmat-tetramer-tratioL0.txt” and “rotmat- 
tetramer-tratioO.5.txt”, respectively. The first, second, 
third columns are the indices i. j. and the entry Uij, re¬ 
spectively. The graphs are numbered from Nl = 6 in 
decreasing order and from No. 1 in increasing order as 
“Graph 1, 2, •••”. 


Appendix B: Average sign computed for the trimer 
and tetramer models 


Figures 11 and [12] show the values of the average sign 
computed for all the graphs of the trimer and tetramer 
models, respectively. See Figs, [l] and [6] for the structures 
of the graphs. 


Appendix C: Distribution function of the expansion 
order in continuous-time QMC 

We consider a quantum system whose Hamiltonian is 
given 


H = H 0 + XHi (A > 0). 


(Cl) 
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FIG. 11. (Color online) Average sign computed for the trimer model: (a) t'/t = 2 and (b) 0.5 at /3 = 50. (Right panel) 
t'/t = 2 and (Left panel) t'/t = 0.5. We show the average sign computed for each basis. 


In a continuous-time QMC algorithm, we expand the par¬ 
tition function as 


Z( A) = Tr[e~P H ] 



xTr[(-l) n e" ( ^" T " )ffo i?i • • • Hie -Tlif °].(C2) 
Practically, we decompose Hi as 


a 


depending on the basis to represent Hi. Then, the par¬ 
tition function is given by 


Z( A) = 


n—0 ‘ 


dri 


dT n \ n w c 


(C4) 


where 


w Cn = Tr[(-1 ) n e-U- Tn ' H °Hi ai ■ ■ ■ Hi an e~ T ' H °}. (C5) 


In a QMC simulation, we sample the partition function 
using w(cn) as MC weight. In the following, we take 
A = 1. The distribution function of n is given by 


P(n) 


£ Cn KJ 

E^oEcJ^cJ’ 


(C3) 


(C6) 
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FIG. 12. (Color online) Average sign computed for the tetramer model at /3 = 50: (Right panel) t'/t = 1 and (Left panel) 
f'/f = 0.5. We show the average sign computed for each basis. 


where the summation is taken over all n-th order 
diagrams. The average sign at each n is defined as fol¬ 
lows: 


(sign) n 


E Crl sign(w> c J|u; c J 

E c „ KJ 


E Cri u, c n 

E c >c„| 


(C7) 


P{n) and (sign)™ generally depe nd on the way of the de¬ 
composition of Hi in Eq. (C3l. However, if (sign)™ is 
independent of n, P{n) does not depend on the decom¬ 
position of Hi. Assuming (sign)™ = C (C is a constant), 


Eq. (C 61 can be written as 


P{n) 


C 1 Ec„ W c n 
c_1 EEoE Cm »c, 

E c „ W Cn 

E OO v—' 

m =0 Ec m 


1 d n Z(\) 
n\ d\ n 


| a=o 


Z{\ = 1) 


(C8) 
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One can easily see that the last line in Eq. (C 81 depends 
only on how we decompose H into Hq and Hi , and hence 
does not depend on the basis in which Hi is expressed. 


Appendix D: Coulomb interaction for bonding and 
anti-bonding orbitals 


Here, we derive the Coulomb interaction for bonding 
and anti-bonding orbitals. We consider a two-site model 
with the onsite Coulomb interaction U whose Hamilto¬ 
nian is given in the site basis as follows: 

2 

H = Uj2nitn ii , (Dl) 

i =1 

where n ia = c\ a c ilT . 

This Hamiltonian can be expressed in the form of a 
Slater-Kanamori interaction as 

1 2 

^ = 2 ^ 1 ^ ', Uq0a'l!l l C < a<T C Ba' C @'< 7 ' Ca '' J ’ (^ 2 ) 

a.[3a.'/3' crcr' 

with the intra-orbital repulsion U aaaa = U, the inter¬ 
orbital repulsion U a / 3 a p = U' = 0, the exchange inter¬ 
action U a pp a = Jh = 0, and the pair hopping Uaapp = 
4 = 0 (a ^/3). 

We now consider the following basis rotation (0 < 9 < 
71 -/ 4 ): 


d a = J2 V ba C b, 

b 

(D3) 

d t = Yl Vba °b > 

b 

(D4) 

v f cos 9 — sin 9 \ 
l sind cos 9 ) ' 

(D5) 


The site basis is given by 9 = 0, while 9 = 7r/4 corre¬ 
sponds to a dimer b asis . 

Substituting Eq. ( |D5| ) into Eq. (D2), one obtains 

'y ' ^ ' Ugfla' 0' c a<T C ])cT ,C P' rr ' Ca ' a 


OLD a.' O' (TO 


0 y ] y ] Uijkld^dj^/dia' dka- 

ijkl oo' 


The Coulomb matrix in the rotated basis reads 


(D6) 


Uijk i = UafJaV'VZVjjVa'kVfn, (D7) 


where i, j, fc, l are the index of orbitals in the rotated 
basis. The nonzero elements in Uijki are 

U = U iiU = (cos 4 9 + sin 4 9)U, (D8) 

U' = Uijij = (2 cos 2 9 sin 2 9)U, (D9) 

J H = Uijji = (2 cos 2 9 sin 2 0)U, (D10) 

Jh = U Ujj = (2 cos 2 9 sin 2 9)U, (Dll) 

where i ^ j. For the dimer basis (9 = 7t/ 4), one obtains 
U = U' = J H = Jh = U/ 2 . 


Appendix E: Improved double vertex update 


In Ref.[32i an efficient CT-INT algorithm, the so-called 
“double vertex update”, has been proposed. The algo¬ 
rithm has been developed to efficiently deal with non¬ 
density type interactions (spin-flip and pair-hopping in¬ 
teractions). In this update, one inserts or removes a pair 
of two vertices corresponding to spin-flip or pair-hopping 
interactions simultaneously. In the present study, we em¬ 
ploy this scheme to deal with the non-density-type inter¬ 
actions generated by the basis transformation. However, 
we found that the acceptance ratio becomes very low at 
large /? if we insert a pair of vertices at two imaginary 
times randomly picked from the whole interval [0, /3) as 
in Ref. 132, This is because the acceptance ratio drops 
exponentially with the time difference between the two 
vertices (see Fig. 13). 

To overcome this problem, we improve the scheme by 
increasing the proposal ratio for double-vertex insertions 
with short time differences. Let us give a detailed expla¬ 
nation by taking the update of two spin-flip interactions 
as an example. The scheme for the pair-hopping inter¬ 
actions can be formulated in an analogous way. For the 
spin-flip interactions, there are two different types of ver¬ 


tices, Jd' a ^dp^d a ±dpf and Jd'g^d/ a ^dpid a f, which appear 
in pairs during the MC simulation. That is, the pertur¬ 
bation orders of these two types of vertices are always 
the same, which we denote by m s . 

In the improved scheme, for the insertion, we pick an 
imaginary time t± randomly in [0, J] for one type of the 
spin-flip interactions. Then, we propose the second imag¬ 
inary time t ' 2 for the other type of the spin-flip interaction 
with a probability proportional to p(Ar) given by 


p(At) = e~ aAr + e - Q ^" Ar) + b (El) 


with At = t ! 2 — r-| |. The positive constants a and b 
should be optimized case by case. For example, the op¬ 
timal value a may depend sensitively on the interaction 
strength. On the other hand, b can generally be chosen 
small, b ~ 10 -3 -10 -2 . p(Ar) is defined in the region 
At £ [0, /?]. One can see that p(Ar) has a larger value 
around At ~ 0 and At ~ J, i.e., we propose the in¬ 
sertion of two vertices with a short time difference more 
frequently. The proposal weight for the insertion of two 
vertices with time difference At is 

-Padd(m s -A m s Tl) = f^) 2 7T ,(E2) 

V P J ^ (1 - exp a P) + b 

where the factor [/L (J — exp _a/3 ) + 6] originates from 
the integral ^ J^ p(r)dT. 

In the removal update, we pick a pair of spin-flip ver¬ 
tices with a probability proportional to p(At) with At 
being the time difference of the two vertices. Then the 
proposal weight for the removal of the pair of the Jth 
spin-flip vertex of type 1 at tj and Jth spin-flip vertex of 
type 2 at t'j is 

Prm(m s + 1 -A m s ) = p(\n - T'j\)/F ms , (E3) 
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At 


FIG. 13. (Color online) Acceptance rate for the insertion of 
a double vertex as a function of the time difference between 
two vertices (At) computed for the tetramer model at U = 8, 
j3 = 50, and t = t' = 1. The data are symmetrized about /3J 2 
and are averaged over spin-flip interactions and pair-hopping 
interactions. 


Then the acceptance rate for the improved double ver¬ 
tex update becomes 


P(c n ->c„ +2 ) = min X J| 


detA a (c n j 2 ) 


det A a {c n ) 


1 (E5) 


for the insertion update and 


P(c„ +2 -»c n ) = min | — 


det A a {cn) 


det Ag. (c n -{_ 2 ) 


, 1 (E6) 


for the removal update. Here, det A determines the 
weight of the configuration in the CT-INT scheme, which 
originates from the Wick’s decomposition of the operator 
series pH n is the number of vertices in the configuration 
including both the density-type and non-density-type in¬ 
teractions. The factor X is given by 


V _ / jj \2 Prm( TO s + l ~^ TO s) 

P&dd(jris m s + l) 


m 2 

^m s +l 


2 

a/3 


(1 - exp _0/3 ) + b 


(E7) 


where F ms is calculated as 

ra s 

F m s = ^2 P(\ r i - T ’j\), (E4) 

i,j—1 

with {Tj} ({rj}) being the imaginary times for the exist¬ 
ing spin-flip vertices of type 1 (type 2). 


For the tetramer, without the improvement presented 
here, the acceptance ratio for the double vertex update 
is as small as about 0.003 for /3 = 50, U = 8, and t = 
t' = 1. We obtained a ~ 5.8 and b ~ 0.002 by fitting the 
accepta nce rate in Fig. 13 using the expression for p(r) 
in Eq. (Ell for 0 < r < 15. We found that with these 
parameters, the acceptance rate is dramatically improved 
from ~ 0.003 to ~ 0.1. 




















